| sample1 | sample2 | sample3 | sample4 | sample5 | sample6 | |
|---|---|---|---|---|---|---|
| gene1 | 85 | 76 | 103 | 107 | 140 | 124 |
| gene2 | 1 | 6 | 11 | 6 | 1 | 4 |
| gene3 | 80 | 98 | 39 | 82 | 97 | 113 |
| gene4 | 92 | 83 | 59 | 85 | 100 | 98 |
| gene5 | 36 | 24 | 18 | 50 | 22 | 15 |
| gene6 | 0 | 0 | 1 | 4 | 2 | 3 |
1 Análisis de expresión génica diferencial (RNA-Seq)
El análisis de la expresión génica diferencial (RNA-Seq) es una de las técnicas bioinformáticas que más auge ha experimentado en la última década, ya que permite cuantificar la expresión de un gen a nivel celular y comparar estas expresiones entre distintos grupos experimentales para identificar asociaciones entre genes y factores experimentales (como enfermedades o tratamientos).
El proceso de un análisis genético diferencial tiene distintas etapas, que van desde la obtención del las células a estudiar, pasando por la secuenciación del RNA que contienen, hasta el análisis estadístico de las frecuencias de expresión de cada gen.


En este tutorial nos centraremos en el análisis estadístico de las frecuencias de expresión de los genes. Para ello partiremos de una matriz con las frecuencias o conteos de expresión de cada gen en el conjunto de muestras analizadas. Las filas de esta matriz corresponden a los genes observados y las columnas a las muestras analizadas.
La mayor parte de los paquetes de Bioconductor que suelen encapsular esta matriz de frecuencias en un objeto del tipo SingleCellExperiment, que además de la matriz de frecuencias de expresión de los genes incluye también los grupos experimentales a los que pertenecen las muestras estudiadas y otra información que se va generando a medida que progresa el análisis.

SingleCellExperimentEl paquete SingleCellExperiment de Bioconductor define esta estructura de datos, y puede instalarse mediante el siguiente código
BiocManager::install('SingleCellExperiment')Existen multitud de paquetes en el repositorio Bioconductor para realizar el análisis de expresión génica diferencial. En este tutorial veremos dos de los más usados: EdgeR y DesSeq2.
1.1 Análisis de expresión génica diferencial con EdgeR
1.1.1 Estructura de datos
El paquete EdgeR utiliza la estructura de datos DGEList para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se utiliza la función
DGEList(frec, group = grupo): Crea una estructura del tipoDGEListcon la matriz de frecuencias de expresión génicafrec(con genes en filas y muestras en columnas) y el vectorgrupocon los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.
library(edgeR)
library(DEFormats)
frec <- simulateRnaSeqData()
grupo <- rep(c("A", "B"), each = 3)
dge <- DGEList(frec, group = grupo)
dgeAn object of class "DGEList"
$counts
sample1 sample2 sample3 sample4 sample5 sample6
gene1 85 76 103 107 140 124
gene2 1 6 11 6 1 4
gene3 80 98 39 82 97 113
gene4 92 83 59 85 100 98
gene5 36 24 18 50 22 15
995 more rows ...
$samples
group lib.size norm.factors
sample1 A 42832 1
sample2 A 40511 1
sample3 A 39299 1
sample4 B 43451 1
sample5 B 40613 1
sample6 B 43662 1
Esta estructura de datos es una lista con dos atributos. El atributo counts contiene la matriz de frecuencias de expresión génica, y el atributo samples es un data frame información sobre las muestras estudiadas. Cada fila de este data frame se corresponde con una columna de la matriz de frecuencias y contiene las siguientes columnas
| Columna | Descripción |
|---|---|
group |
Grupo experimental al que pertenece la muestra. |
lib.size |
tamaño de la librería (suma de frecuencias de la muestra). |
norm.factors |
Factor de normalización. |
La estructura de datos DGEList puede contener opcionalmente el atributo genes con anotaciones de los genes observados en las filas de la matriz de frecuencias.
Para convertir esta estructura de datos en una del tipo DESeqDataSet se puede utilizar la función as.DESeqDataSet del paquete DEFormats.
dds <- as.DESeqDataSet(dge)
ddsclass: DESeqDataSet
dim: 1000 6
metadata(1): version
assays(1): counts
rownames(1000): gene1 gene2 ... gene999 gene1000
rowData names(0):
colnames(6): sample1 sample2 ... sample5 sample6
colData names(3): group lib.size norm.factors
Ejemplo 1.1 Veamos cómo crear la estructura de datos correspondiente al estudio de Sheridan JM, et al. (2015) donde se obtuvieron datos de tres poblaciones de células (basal, progenitor luminal (LP) y luminal maduro (ML)) seleccionadas de las glándulas mamarias de ratones vírgenes hembra, cada una por triplicado. Los datos pueden obtenerse del repositorio Gene Expression Omnibus (GEO) mediante el identificador GSE63310.
En primer lugar descargamos los ficheros con las librerías de frecuencias.
# Descarga de datos
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="datos/GSE63310_RAW.tar", mode="wb")
utils::untar("datos/GSE63310_RAW.tar", exdir = "datos")
files <- paste0("datos/", c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
"GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt"))
#for(i in paste0("datos/", files, ".gz"))
# R.utils::gunzip(i, overwrite=TRUE)A continuación creamos la estructura de datos DGEList.
# Creación de la estructura de datos DGEList
dge <- readDGE(files, columns=c(1,3))
dgeAn object of class "DGEList"
$samples
files group lib.size
datos/GSM1545535_10_6_5_11 datos/GSM1545535_10_6_5_11.txt 1 32863052
datos/GSM1545536_9_6_5_11 datos/GSM1545536_9_6_5_11.txt 1 35335491
datos/GSM1545538_purep53 datos/GSM1545538_purep53.txt 1 57160817
datos/GSM1545539_JMS8-2 datos/GSM1545539_JMS8-2.txt 1 51368625
datos/GSM1545540_JMS8-3 datos/GSM1545540_JMS8-3.txt 1 75795034
datos/GSM1545541_JMS8-4 datos/GSM1545541_JMS8-4.txt 1 60517657
datos/GSM1545542_JMS8-5 datos/GSM1545542_JMS8-5.txt 1 55086324
datos/GSM1545544_JMS9-P7c datos/GSM1545544_JMS9-P7c.txt 1 21311068
datos/GSM1545545_JMS9-P8c datos/GSM1545545_JMS9-P8c.txt 1 19958838
norm.factors
datos/GSM1545535_10_6_5_11 1
datos/GSM1545536_9_6_5_11 1
datos/GSM1545538_purep53 1
datos/GSM1545539_JMS8-2 1
datos/GSM1545540_JMS8-3 1
datos/GSM1545541_JMS8-4 1
datos/GSM1545542_JMS8-5 1
datos/GSM1545544_JMS9-P7c 1
datos/GSM1545545_JMS9-P8c 1
$counts
Samples
Tags datos/GSM1545535_10_6_5_11 datos/GSM1545536_9_6_5_11
497097 1 2
100503874 0 0
100038431 0 0
19888 0 1
20671 1 1
Samples
Tags datos/GSM1545538_purep53 datos/GSM1545539_JMS8-2
497097 342 526
100503874 5 6
100038431 0 0
19888 0 0
20671 76 40
Samples
Tags datos/GSM1545540_JMS8-3 datos/GSM1545541_JMS8-4
497097 3 3
100503874 0 0
100038431 0 0
19888 17 2
20671 33 14
Samples
Tags datos/GSM1545542_JMS8-5 datos/GSM1545544_JMS9-P7c
497097 535 2
100503874 5 0
100038431 1 0
19888 0 1
20671 98 18
Samples
Tags datos/GSM1545545_JMS9-P8c
497097 0
100503874 0
100038431 0
19888 0
20671 8
27174 more rows ...
Añadimos la información del diseño experimental, en este caso el grupo experimental y el lote, al data frame de las muestras.
colnames(dge) <- substring(colnames(dge), 18, nchar(colnames(dge)))
# Añadimos el grupo experimental
dge$samples$group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
# Añadimos el lote
dge$samples$lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
dgeAn object of class "DGEList"
$samples
files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt LP 32863052 1 L004
9_6_5_11 datos/GSM1545536_9_6_5_11.txt ML 35335491 1 L004
purep53 datos/GSM1545538_purep53.txt Basal 57160817 1 L004
JMS8-2 datos/GSM1545539_JMS8-2.txt Basal 51368625 1 L006
JMS8-3 datos/GSM1545540_JMS8-3.txt ML 75795034 1 L006
JMS8-4 datos/GSM1545541_JMS8-4.txt LP 60517657 1 L006
JMS8-5 datos/GSM1545542_JMS8-5.txt Basal 55086324 1 L006
JMS9-P7c datos/GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
JMS9-P8c datos/GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
$counts
Samples
Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
497097 1 2 342 526 3 3 535 2
100503874 0 0 5 6 0 0 5 0
100038431 0 0 0 0 0 0 1 0
19888 0 1 0 0 17 2 0 1
20671 1 1 76 40 33 14 98 18
Samples
Tags JMS9-P8c
497097 0
100503874 0
100038431 0
19888 0
20671 8
27174 more rows ...
A continuación anotamos los genes de las filas de la matriz de frecuencias. Para ello debe utilizarse un paquete específico para el genoma de la especie de donde proviene el RNA (Mus.muculus en este caso para el genoma del ratón).
library(Mus.musculus)
geneid <- rownames(dge)
genes <- select(Mus.musculus, keys = geneid, columns = c("SYMBOL", "TXCHROM"), keytype = "ENTREZID")
# Eliminamos duplicidad de algunos genes
genes <- genes[!duplicated(genes$ENTREZID),]
dge$genes <- genes
dgeAn object of class "DGEList"
$samples
files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt LP 32863052 1 L004
9_6_5_11 datos/GSM1545536_9_6_5_11.txt ML 35335491 1 L004
purep53 datos/GSM1545538_purep53.txt Basal 57160817 1 L004
JMS8-2 datos/GSM1545539_JMS8-2.txt Basal 51368625 1 L006
JMS8-3 datos/GSM1545540_JMS8-3.txt ML 75795034 1 L006
JMS8-4 datos/GSM1545541_JMS8-4.txt LP 60517657 1 L006
JMS8-5 datos/GSM1545542_JMS8-5.txt Basal 55086324 1 L006
JMS9-P7c datos/GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
JMS9-P8c datos/GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
$counts
Samples
Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
497097 1 2 342 526 3 3 535 2
100503874 0 0 5 6 0 0 5 0
100038431 0 0 0 0 0 0 1 0
19888 0 1 0 0 17 2 0 1
20671 1 1 76 40 33 14 98 18
Samples
Tags JMS9-P8c
497097 0
100503874 0
100038431 0
19888 0
20671 8
27174 more rows ...
$genes
ENTREZID SYMBOL TXCHROM
1 497097 Xkr4 chr1
2 100503874 Gm19938 <NA>
3 100038431 Gm10568 <NA>
4 19888 Rp1 chr1
5 20671 Sox17 chr1
27174 more rows ...
1.1.2 Preprocesamiento
Una vez preparada la estructura de datos, el siguiente paso es el preprocesamiento de datos. Normalmente comprende las siguientes tareas:
- Normalización de las frecuencias.
- Eliminación de los genes con poca expresión.
- Normalización de las distribuciones de frecuencias de expresión génicas.
1.1.2.1 Normalización de las frecuencias
El objetivo principal es normalizar las frecuencias de expresión génica para eliminar el efecto sobre las frecuencias de factores como la profundidad de secuenciado (a mayor profundidad de secuenciado mayores frecuencias) o el tamaño de los genes (a mayor tamaño mayores frecuencias). Generalmente se usan las siguientes transformaciones
Frecuencias por millón (CPM). Se divide cada frecuencia por el tamaño en millones de la librería de la muestra a la que pertenece. Se realiza con la función
cpm(dge).Logaritmo en base 2 de las frecuencias por millón (log2-CPM). Se calcula a partir de la anterior mediante la fórmula \(\log_2(CPM+\frac{2}{L})\), donde \(L\) es la longitud media de las librerías en millones. El término \(\frac{2}{L}\) que se añade permite calcular el logaritmo para frecuencias nulas (ya que el logaritmo de 0 no existe). Se realiza con la función
cmp(dge, log = TRUE).Lecturas por kilobase de transcripción por millón (RPKM): Se realiza con la función
rpkm(dge, longitud_genes)pasando la longitud de los genes.Fragmentos por kilobase de transcripción por millón (FPKM).
Las dos primeras no tienen en cuenta el tamaño de los genes, pero suelen usarse para la expresión génica diferencial donde el tamaño de los genes se supone constante en las distintas muestras.
Ejemplo 1.2 Siguiendo con el ejemplo anterior, calculamos las frecuencias por millón y el logaritmo de las frecuencias por millón.
cpm <- cpm(dge)
summary(cpm) 10_6_5_11 9_6_5_11 purep53 JMS8-2
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.578 Median : 0.736 Median : 0.892 Median : 0.895
Mean : 36.793 Mean : 36.793 Mean : 36.793 Mean : 36.793
3rd Qu.: 19.536 3rd Qu.: 23.546 3rd Qu.: 24.221 3rd Qu.: 23.341
Max. :27807.947 Max. :11546.719 Max. :7951.408 Max. :7389.433
JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.699 Median : 0.711 Median : 0.908 Median : 0.845
Mean : 36.793 Mean : 36.793 Mean : 36.793 Mean : 36.793
3rd Qu.: 23.827 3rd Qu.: 19.928 3rd Qu.: 21.439 3rd Qu.: 24.260
Max. :7955.680 Max. :29572.361 Max. :9376.973 Max. :7865.350
JMS9-P8c
Min. : 0.000
1st Qu.: 0.000
Median : 0.752
Mean : 36.793
3rd Qu.: 21.594
Max. :16500.710
lcpm <- cpm(dge, log = TRUE)
summary(lcpm) 10_6_5_11 9_6_5_11 purep53 JMS8-2
Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
Median :-0.6847 Median :-0.3589 Median :-0.09513 Median :-0.0901
Mean : 0.1714 Mean : 0.3312 Mean : 0.43559 Mean : 0.4089
3rd Qu.: 4.2913 3rd Qu.: 4.5601 3rd Qu.: 4.60081 3rd Qu.: 4.5475
Max. :14.7632 Max. :13.4952 Max. :12.95700 Max. :12.8513
JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
Median :-0.4281 Median :-0.4064 Median :-0.07152 Median :-0.1704
Mean : 0.3225 Mean : 0.2529 Mean : 0.40428 Mean : 0.3708
3rd Qu.: 4.5772 3rd Qu.: 4.3199 3rd Qu.: 4.42513 3rd Qu.: 4.6031
Max. :12.9578 Max. :14.8520 Max. :13.19491 Max. :12.9413
JMS9-P8c
Min. :-4.5074
1st Qu.:-4.5074
Median :-0.3300
Mean : 0.2749
3rd Qu.: 4.4355
Max. :14.0102
1.1.2.2 Eliminación de genes con poca expresión
Aunque el objetivo del análisis de la expresión génica diferencial es detectar los genes que se expresan en un grupo experimental en comparación con los otros, cuando un gen no se expresa en ninguna de las muestras no tiene interés para el estudio y puede eliminarse.
Existen diferentes criterios para eliminar los genes con poca expresión. El paquete EdgeR incorpora la siguiente función
filterByExpr(dge): Realiza un filtro de la estructura de datosdgelos genes con poca expresión. Por defecto devuelveTRUEpara cualquier gen con una frecuencia mayor o igual que 10 en al menos un número de muestras igual al del menor grupo experimental.
Ejemplo 1.3 Veamos cuántos genes no tienen expresión en ninguna muestra del ejemplo anterior.
# Número de genes con expresión nula.
table(rowSums(dge$counts) == 0)
FALSE TRUE
22026 5153
El siguiente gráfico muestra la distribución del logaritmo en base 2 de las frecuencias por millón.
# Definimos una función para dibujar la distribución del logaritmo de las frecuencias por millón.
dist_logcpm <- function(lcpm) {
lcpm |>
as.tibble() |>
pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
ggplot(aes(x = LogCPM, color = Muestra)) +
geom_density() +
labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}
dist_logcpm(lcpm)Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Como se observa ha una porcentaje significativo de genes que tiene una expresión muy baja (valores negativos).
A continuación eliminamos de la estructura de datos los genes con poca expresión.
filtro <- filterByExpr(dge)
dge <- dge[filtro, , keep.lib.sizes = F]
dgeAn object of class "DGEList"
$samples
files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt LP 32857304 1 L004
9_6_5_11 datos/GSM1545536_9_6_5_11.txt ML 35328624 1 L004
purep53 datos/GSM1545538_purep53.txt Basal 57147943 1 L004
JMS8-2 datos/GSM1545539_JMS8-2.txt Basal 51356800 1 L006
JMS8-3 datos/GSM1545540_JMS8-3.txt ML 75782871 1 L006
JMS8-4 datos/GSM1545541_JMS8-4.txt LP 60506774 1 L006
JMS8-5 datos/GSM1545542_JMS8-5.txt Basal 55073018 1 L006
JMS9-P7c datos/GSM1545544_JMS9-P7c.txt ML 21305254 1 L008
JMS9-P8c datos/GSM1545545_JMS9-P8c.txt LP 19955335 1 L008
$counts
Samples
Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
497097 1 2 342 526 3 3 535 2
20671 1 1 76 40 33 14 98 18
27395 431 771 1368 1268 1564 769 818 468
18777 768 1722 2517 1923 3865 1888 1830 1246
21399 810 977 2472 1870 2251 1716 1932 756
Samples
Tags JMS9-P8c
497097 0
20671 8
27395 342
18777 693
21399 619
16619 more rows ...
$genes
ENTREZID SYMBOL TXCHROM
1 497097 Xkr4 chr1
5 20671 Sox17 chr1
6 27395 Mrpl15 chr1
7 18777 Lypla1 chr1
9 21399 Tcea1 chr1
16619 more rows ...
lcpm <- cpm(dge, log = TRUE)
dist_logcpm(lcpm)
1.1.2.3 Normalización de las distribuciones de frecuencias de expresión génicas
Durante el proceso de secuenciación puede haber factores externos no biológicos que afecten a la expresión de los genes. Por ejemplo, la muestras del primer lote pueden tener mayores frecuencias que las del segundo lote. Como se supone que las muestras replicadas deben tener un rango y una distribución de frecuencias similares, en esta etapa se aplica otro procedimiento de normalización para garantizar que la distribución de frecuencias de cada muestra es similar. Para ello el paquete edgeR utiliza el método de las medias recortadas de los valores M por medio de la función
calcNormFactors(dge, method = "TMM"): Calcular los factores de escalado del tamaño de las librerías de cada muestra. Estos factores se guardan automáticamente en el data frame con la información de las muestrasdge$samples.
Ejemplo 1.4 Siguiendo con el mismo ejemplo, calculamos los factores de escalado para cada muestra.
dge <- calcNormFactors(dge, method = "TMM")
dge$samples$norm.factors[1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959
[8] 1.0861026 0.9839203
A continuación se muestran los diagramas de cajas de las distribuciones normalizadas tras aplicar los factores de escalado.
box_logcpm <- function(lcpm){
lcpm |>
as.tibble() |>
pivot_longer(everything(), names_to = "Muestra", values_to = "LogCPM") |>
ggplot(aes(x = Muestra, y = LogCPM, fill = Muestra)) +
geom_boxplot() +
labs(title = "Distribución del logaritmo en base 2 de las frecuencias por millón")
}
box_logcpm(lcpm)
Como se puede apreciar, después de la normalización, todas las muestras presentan una distribución de frecuencias similar.
1.1.3 Análisis exploratorio
Una vez preprocesados los datos comienza el análiss estadístico propiamente dicho. En una primera fase se suele realizar un análisis exploratorio de los datos que suele incluir los siguientes análisis:
- Escalado multidimensional (análisis de componentes principales).
1.1.3.1 Escalado multidimensional
El escalado multidimensional mediante componente principales permite ver qué muestras son similares en cuando a la distribución de frecuencias de expresión génica. Los componentes principales son combinaciones lineales de los genes de la matriz de frecuencias que son ortogonales entre sí. El primer componente principal recoge la dimensión con mayor variabilidad de las frecuencias. El segundo recoge la dimensión don la mayor variabilidad no explicada por el primer componente principal y así sucesivamente. Normalmente los dos primeros componentes principales suelen recoger un porcentaje bastante alto de la variabilidad de las frecuencias. Al representar las muestras en estas dos dimensiones, las muestras más próximas entre sí, presentan una distribución de frecuencias de expresión génica similar. Para realizar esta representación se puede utilizar la siguiente función del paquete limma:
plotMDS(dge): Realiza un diagrama de componentes principales de la matriz de frecuencias de la estructura de datosdge`.
Alternativamente, se pueden calcular los componentes principales mediante la función prcomp del paquete stats y luego usar la función autoplot del paquete ggfortify para dibujar el diagrama de componentes principales.
Ejemplo 1.5 A continuación se muestra cómo obtener el diagrama de componentes principales de la matriz de frecuencias de nuestro ejemplo.
plotMDS(dge, col = as.numeric(dge$samples$group), main = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")
library(ggfortify)
autoplot(prcomp(t(lcpm)), data = dge$samples, color = "group", loadings = TRUE, loadings.label = TRUE) +
labs(title = "Diagrama de componentes principales del logaritmo en base 2 de las frecuencias por millón")
Como se puede apreciar las muestras de cada grupo experimental aparecen agrupadas. La mayor diferencia (a lo largo del primer componente principal) se da entre el grupo basal y los otros dos grupos. Esto indica que cuando se realice el análisis de expresión génica diferencial habrá bastantes genes con diferencias de expresión significativa entre el grupo basal y los otros dos, mientras que no habrá tantos al comparar los grupos LP y ML.
Otra opción muy interesante es el paquete Glimma que permite dibujar un diagrama de componentes principales interactivo mediante la función
glMDSPlot(lcpm): Dibuja un diagrama de componentes principales interactivo de la matriz de frecuenciaslcpm.
library(Glimma)
#glMDSPlot(lcpm, groups = dge$samples[,c(2,5)])
dds <- DESeqDataSetFromMatrix(
countData = dge$counts,
colData = dge$samples,
rowData = dge$genes,
design = ~group
)converting counts to integer mode
glimmaMDS(dds)1.1.4 Análisis de expresión génica diferencial
La última etapa del análisis es la detección de los genes con una expresión significativamente diferente entre los grupos experimentales. Para ello se suelen utilizar modelos lineales o modelos lineales generalizados. En general, la estimación de los parámetros del modelo ajustado depende de la distribución teórica usada para modelizar las frecuencias de expresión génica. El paquete limma, por ejemplo, realiza un ajuste de modelo lineal suponiendo que las distribuciones de las variables implicadas son normales, mientras que el paquete edgeR modeliza las frecuencias de expresión de los genes observadas mediante la distribución binomial negativa.
El primer paso es definir la matriz del diseño del experimento, que define los grupos experimentales a comparar.
diseño <- model.matrix(~ 0 + dge$samples$group + dge$samples$lane)
colnames(diseño) <- gsub("dge\\$samples\\$group", "", colnames(diseño))
colnames(diseño) <- gsub("dge\\$samples\\$", "", colnames(diseño))
diseño Basal LP ML laneL006 laneL008
1 0 1 0 0 0
2 0 0 1 0 0
3 1 0 0 0 0
4 1 0 0 1 0
5 0 0 1 1 0
6 0 1 0 1 0
7 1 0 0 1 0
8 0 0 1 0 1
9 0 1 0 0 1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"
attr(,"contrasts")$`dge$samples$lane`
[1] "contr.treatment"
Para definir los contrastes de ajuste por pares, el paquete limma utiliza la función
makeContrast(pares)
contr.matrix <- makeContrasts(
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(diseño))
contr.matrix Contrasts
Levels BasalvsLP BasalvsML LPvsML
Basal 1 1 0
LP -1 0 1
ML 0 -1 -1
laneL006 0 0 0
laneL008 0 0 0
En general en las distribuciones de frecuencias de expresión génica, se ha observado que la varianza no es independiente de la media. Los métodos que modelizan las frecuencias mediante el modelo binomial negativo asumen una relación cuadrática entre la media y la varianza. Con el paquete limma el ajuste lineal se realiza sobre el logaritmo en base 2 de las frecuencias por millón (log-CPM) que se suponen tienen una distribución normal. La relación entre la media y la varianza se realiza automáticamente mediante la función voom.
par(mfrow=c(1,2))
v <- voom(dge, diseño, plot=TRUE)
vAn object of class "EList"
$genes
ENTREZID SYMBOL TXCHROM
1 497097 Xkr4 chr1
5 20671 Sox17 chr1
6 27395 Mrpl15 chr1
7 18777 Lypla1 chr1
9 21399 Tcea1 chr1
16619 more rows ...
$targets
files group lib.size norm.factors lane
10_6_5_11 datos/GSM1545535_10_6_5_11.txt LP 29387429 0.8943956 L004
9_6_5_11 datos/GSM1545536_9_6_5_11.txt ML 36212498 1.0250186 L004
purep53 datos/GSM1545538_purep53.txt Basal 59771061 1.0459005 L004
JMS8-2 datos/GSM1545539_JMS8-2.txt Basal 53711278 1.0458455 L006
JMS8-3 datos/GSM1545540_JMS8-3.txt ML 77015912 1.0162707 L006
JMS8-4 datos/GSM1545541_JMS8-4.txt LP 55769890 0.9217132 L006
JMS8-5 datos/GSM1545542_JMS8-5.txt Basal 54863512 0.9961959 L006
JMS9-P7c datos/GSM1545544_JMS9-P7c.txt ML 23139691 1.0861026 L008
JMS9-P8c datos/GSM1545545_JMS9-P8c.txt LP 19634459 0.9839203 L008
$E
Samples
Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5
497097 -4.292165 -3.856488 2.5185849 3.2931366 -4.459730 -3.994060 3.2869677
20671 -4.292165 -4.593453 0.3560126 -0.4073032 -1.200995 -1.943434 0.8442767
27395 3.876089 4.413107 4.5170045 4.5617546 4.344401 3.786363 3.8990635
18777 4.708774 5.571872 5.3964008 5.1623650 5.649355 5.081611 5.0602470
21399 4.785541 4.754537 5.3703795 5.1220551 4.869586 4.943840 5.1384776
Samples
Tags JMS9-P7c JMS9-P8c
497097 -3.2103696 -5.295316
20671 -0.3228444 -1.207853
27395 4.3396075 4.124644
18777 5.7513694 5.142436
21399 5.0308985 4.979644
16619 more rows ...
$weights
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1.079413 1.332986 19.826915 20.273253 1.993686 1.395853 20.494977
[2,] 1.170357 1.456380 4.804866 8.660025 3.612508 2.626870 8.760149
[3,] 20.219073 25.573792 30.434759 28.528310 31.352260 25.743247 28.722497
[4,] 26.947557 32.505933 33.583128 33.232125 34.231754 32.354158 33.334340
[5,] 26.610864 28.501638 33.645479 33.206374 33.573492 31.996623 33.308490
[,8] [,9]
[1,] 1.107780 1.079413
[2,] 3.211473 2.541942
[3,] 21.200072 16.657930
[4,] 30.348630 24.259801
[5,] 25.171513 23.573305
16619 more rows ...
$design
Basal LP ML laneL006 laneL008
1 0 1 0 0 0
2 0 0 1 0 0
3 1 0 0 0 0
4 1 0 0 1 0
5 0 0 1 1 0
6 0 1 0 1 0
7 1 0 0 1 0
8 0 0 1 0 1
9 0 1 0 0 1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`dge$samples$group`
[1] "contr.treatment"
attr(,"contrasts")$`dge$samples$lane`
[1] "contr.treatment"
vfit <- lmFit(v, diseño)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
A partir del modelo estimado se pueden obtener los genes subexpresados y sobreexpresados mediante el p-valor ajustado, tomando un nivel de significación del 5%.
summary(decideTests(efit)) BasalvsLP BasalvsML LPvsML
Down 4648 4927 3135
NotSig 7113 7026 10972
Up 4863 4671 2517
Para la comparación entre los niveles de expresión en basal y LP, se encontró que 4648 genes están regulados negativamente en basal en relación con LP y 4.863 genes están regulados al alza en basal en relación con LP, un total de 9.511 genes significativamente diferenciados. Por otro lado, hay un total de 9598 genes significativamente diferenciados entre basal y ML (4927 genes regulados negativamente y 4671 regulados positivamente), y un total de 5652 genes significativamente diferenciados entre LP y ML (3135 regulados negativamente y 2517 regulados positivamente).
En algunos análisis se puede ser más conservador utilizando el log-fold-change (log-FC).
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt) BasalvsLP BasalvsML LPvsML
Down 1632 1777 224
NotSig 12976 12790 16210
Up 2016 2057 190
Para terminar se pueden mostrar los genes con diferencias más significativas.
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)
head(basal.vs.lp) ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value
12759 12759 Clu chr14 -5.455444 8.856581 -33.55508 1.723731e-10
53624 53624 Cldn7 chr11 -5.527356 6.295437 -31.97515 2.576972e-10
242505 242505 Rasef chr4 -5.935249 5.118259 -31.33407 3.081544e-10
67451 67451 Pkp2 chr16 -5.738665 4.419170 -29.85616 4.575686e-10
228543 228543 Rhov chr2 -6.264208 5.485314 -29.07484 5.782844e-10
70350 70350 Basp1 chr15 -6.084738 5.247023 -28.26649 7.267694e-10
adj.P.Val
12759 1.707586e-06
53624 1.707586e-06
242505 1.707586e-06
67451 1.739242e-06
228543 1.739242e-06
70350 1.739242e-06
basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)
head(basal.vs.ml) ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value
242505 242505 Rasef chr4 -6.534668 5.118259 -35.08270 1.226399e-10
53624 53624 Cldn7 chr11 -5.495888 6.295437 -31.68918 2.774140e-10
12521 12521 Cd82 chr2 -4.690834 7.069637 -31.43063 2.913572e-10
20661 20661 Sort1 chr3 -4.931660 6.704459 -30.70113 3.558720e-10
71740 71740 Nectin4 chr1 -5.581283 5.164967 -30.59775 3.718157e-10
12759 12759 Clu chr14 -4.686826 8.856581 -27.95760 7.687544e-10
adj.P.Val
242505 1.236213e-06
53624 1.236213e-06
12521 1.236213e-06
20661 1.236213e-06
71740 1.236213e-06
12759 1.480434e-06
glimmaMA(tfit, coef=1, status=dt, main=colnames(tfit)[1],
side.main="ENTREZID", counts=lcpm, groups=dge$samples$group)External counts supplied using counts argument will be transformed to log-cpm by default. Specify transform.counts='none' to override transformation.
Warning in buildXYData(table, status, main, display.columns, anno, counts, :
count transform requested but not all count values are integers.
El mismo ajuste se puede realizar con el paquete edgeR.
dge <- estimateDisp(dge, diseño)
gfit <- glmFit(dge, diseño)
glrt <- glmLRT(gfit, diseño, contrast = contr.matrix)#glimmaMA(glrt, dge = dge)Y el volcano plot
glimmaVolcano(efit, dge = dge)1.2 Análisis de expresión génica diferencial con DesSeq2
El paquete DESeq2 utiliza la estructura de datos DESeqDataSet para guardar la matriz de frecuencias de expresión génica. Para crear esta estructura se pueden utilizar las siguientes funciones
DESeqDataSetFromMatrix(countData = frec, colData = grupo): Crea una estructura del tipoDESeqDataSetcon la matriz de frecuencias de expresión génicafrec(con genes en filas y muestras en columnas), y el data framegrupocon los grupos experimentales a los que pertenecen las muestras analizadas en las columnas de la matriz de frecuencias.
#library(DESeq2)
#dds <- DESeqDataSetFromMatrix(countData = frec, colData = dge$samples$group, design = ~ group)
#dds1.3 Referencias
Orchestrating Single-Cell Analysis with Bioconductor. (s. f.). Recuperado 12 de junio de 2023, de https://github.com/LTLA/OSCA
Law, Charity, et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR.